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C/^ ■ Abstract 

0^ ■ In many epidemiological models a nonlinear transmission function is used in 

,— i, the form of power law relationship. It is constantly argued that such form reflects 

I population heterogeneities including differences in the mixing pattern, susceptibility, 

Ph ■ and spatial patchiness, although the function itself is considered phenomenological. 

O . Comparison with large-scale simulations show that models with this transmission 

[ function accurately approximate data from highly heterogeneous sources. In this 

qh' note we provide a mechanistic derivation of the power law transmission function, 
I— starting with a simple heterogeneous susceptibles-infectives (SI) model, which is 

I based on a standard mass action assumption. We also consider the simplest SI 

J> ' model with separable mixing and compare our results with known results from the 

9^ ■ literature. 
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X ; 1 Introduction 

It is customary to consider transmission function T{S,I), wliicli describes tlie incidence 
rate, i.e., tlie number of new cases per time unit, as a main component of any epidemio- 
logical model 0, 20]. Here we use usual notations for susceptible and infective individuals 



denoting them as S and / respectively. Assuming that there is no influx of susceptible 
hosts in out model, we can write that 

|sM = -r(s./). 
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Historically the earliest form of the transmission function was a simple bilinear form, 
i.e., T(S,I) = PSI which follows from the assumptions of random contacts, host 
homogeneity, and application of the law of mass action, thereby implying that the contact 
rate of any individual is a linear function of the population size (see [5)] for more details). 
Here (3 > is the transmission coefficient. Under the proportional mixing assumption (the 
contact rate is fixed), the transmission function takes the form T{S, I) = f3SI/N, iV is a 
population size. If the model includes an assumption of the constant population size these 
two transmission functions are virtually the same from any practical viewpoint, whereas 
variable population size can yield dramatically different behaviors (e.g., [2, Isl, 122] ) . 

It was early acknowledged that other than bilinear or proportional mixing transmission 
functions should be used in epidemiological models to provide better fit of the model solu- 
tions to empirical data (see 2d| for a general account of different models for transmission 



functions). 

One of the most widely used functions has the following form: 

TiS,I)=f3SPI'^, p, g>0. (1) 
We will term this transmission function as power law relationship. It was first used in (29l . 



30j in the form T{S, I) = (3S^I "to investigate the consequences of various assumptions 
when the laws are not known" . Severo [25,] considered general form ([T]) where both p and 
q are not equal to one, though he did not give a detailed analysis of the model. Liu et 



al. 1 181 . Il9[ gave a thorough analysis of different compartmental epidemiological models 
with ([T]) and showed that incorporating power law transmission function yields various 
dynamical behaviors not observable in models with bilinear incidence rate, e.g., limit 
cycles and multiple equihbrium points. Additional analysis and details of such models 
can be found in fol. idk. 



Since the first use of the power law transmission function its form was explained on a 
basis of "intrinsic heterogeneity in mixing pattern" of a population under question. The 



exponents p and q were dubbed as "heterogeneity parameters" |25| . but the model itself 
is considered phenomenological and lacking mechanical derivation 20] in contract to, e.g., 
bilinear relationship, which is based on a dubious but well established law of mass action 

i. 

The link between phenomenological power law incidence rate and population hetero- 
geneity was made explicit when it was shown that such mean-field models can provide 
an accurate approximation to network based simulations that include variation in the 



strength, duration, and number of contacts per person. In [26j the transmission function 
was used in the form T{S, I) = jSS^I, whereas full non-linear transmission function was 
implemented in ^]. In both cases it was shown that power law relationship improves 
the accuracy of mean-field model predictions when compared with models with bilin- 
ear transmission function (see also [1] for a review on comparison of homogeneous and 
heterogeneous models). 

In this note we show that power law transmission function can be not only postulated 
but also derived, using a simple heterogeneous SI model. The paper organized as follows. 
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In the next section we formulate a mechanistic heterogeneous SI model from the first 
principles. Section 3 gives a brief exposition of necessary analytical tools. In Section 4 we 
present the main results of the study showing that our heterogeneous model is equivalent 
to a homogeneous one, but with a non-linear transmission function. 



2 Model formulation 



Heterogeneity profoundly affects the dynamics of infection. Differences in contact rates, 
spatial distributions of susceptible hosts, infectiousness and susceptibility of individuals 
have a direct effect on disease dynamics. Here, we specifically look into heterogeneity 
in disease parameters (such as susceptibility) do not touching an important topics of 
heterogeneity mediated by a structured variable, such as explicit space or age structure. 
Our approach is close to the one given in, e.g., p, 0,123] (see also 23] for more details). 



Model I. We start with a generic assumption that the subpopulation of susceptible hosts 
is heterogeneous, and denote s{t,u) the density of susceptibles at time t having parameter 
value uj, which determines susceptibility to a particular disease and varies from individual 
to individual. The total size of the susceptibles is given by S{t) = J^s{t,uj) duj, where 
Q is the set of parameter values. Assuming that the subpopulation of the infectives is 
homogeneous (later we relax this assumption), the contact process is described with the 
law of mass action, and the rate of change in the susceptibles is determined by transmission 
parameter, which is a function of u, we obtain that 

^^s{t,uj) = -(3{uj)s{t,uj)I{t). (2) 

Here I3{uj) incorporates information on the contact rate and the probability of a successful 
contact. 

The change in the infective class is given by 

^I{t) = lit) [ s{t, u) du = mS{t)I{t), (3) 



dt 



n 



where we denote 



f3{t) I ps{t,uj)duj, Ps{t,uj) ^'^^'^^ 



s{t) ■ 

Therefore, (3{t) is the mean value of the function (3{uj) with respect to probability density 
function ps{t,Lj) for any time t. We need the initial conditions for the model (ED, (El): 

s(0, u) = soiuj) = SoPsiO, uj), /(O) = Jo. (4) 

Here Sq, Iq are given numbers, and ps{0,uj) is a given initial distribution of the suscepti- 
bility in the population. 
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We note that formally, after integrating equation with respect to u, we obtain 
a homogeneous SI model with non-constant transmission parameter P{t) which, in its 
turn, depends on the current distribution of susceptibility in the population. If f3(t) is 
known then the problem is solved. Interesting to remark that ad hoc approach to use 
time-dependent transmission coefficient in an SIR model was used to approximate a 
heterogeneous epidemics with a mean- field model jl6| . 



Model II. Let us assume now that not only the susceptibles are heterogeneous for 
some trait that infiuences the disease evolution, but also the infectives are heteroge- 
neous, and consider the simplest possible SI model. Let s(t,a;i) and i{t,uj2) be the 
densities of the susceptibles and infectives respectively, here we assume that the traits 
of the two classes are independent, i.e., /3(ci;i, 0^2) = /?i(co'i)/?2('^2)- The number of sus- 
ceptibles with the trait value ui infected by individuals with trait value UJ2 is given by 
/?i(ct;i)s(t, Ci;i)/32 (^2)^(^1 ^^2); and the total change in the infective class with trait value 
UJ2 is P2{^2)i(t, LJ2) /j^^ /?i(co'i)s(t, Ui) dui, an analogous expression applies to the change in 
the susceptible population. We emphasize that nothing else except for the standard law of 
mass action is supposed to formulate the terms for the change in susceptible and infective 
subpopulations. Combining the above assumptions we obtain the following model: 



d f 

— s{t,UJi) = -Pi{uJi)s{t,UJi) / j32{uJ2)i{t,^2)dU2 

'^t Jn2 

= -Pl{uJ,)s{t,UJ,)P2{t)I{t), 
d f 

—z{t,UJ2) = P2Mi{t,UJ2) J (3i{uJi)s{t,Ui)duJi 
= f32{uJ2)t{t, 002) Mt)S{t). 

Model ([5]) is supplemented with the initial conditions s(0,c<Ji) = SoPs{0,uJi), i{0,u!2) = 

IoPi{0,UJ2). 

In (EJ it is assumed that if an individual having trait value ui was infected by an 
individual with trait value 002 he or she becomes an infective with trait value 002- This 
is a restrictive assumptions which is necessary to apply the main theorem from the next 
section. 

The global dynamics of ([5]), as well as of (l2])-(jll), is simple and is similar to the simplest 
homogeneous SI model. 



Model III. Above we were talking about heterogeneity of the hosts: whether all suscep- 
tible individuals are of the same type with equal susceptibility, and whether all infectious 
individuals have equal ability to infect others. Another aspect of heterogeneity is the pos- 
sible heterogeneous social contact network It is difficult to apply the general theory 
of heterogeneous populations (see below) to such models, however, there is a simple case, 
for which some results can be obtained. 
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Let us assume that n(t,u!) denotes the density of individuals in the population, which 
are making u contacts on average. Every individual can be contacted by another individ- 
ual, which differs is an average number of contact per individual. This situation is usually 
termed as separable mixing. If we denote r the probability of transmission the disease 
given a contact, then, the simplest SI- model with separable mixing can be described by 
the following system: 

d , . , , Luji{t,uj)duj 

^At,w) = -njs{t,u;)^ — — , 

dt J^ujno{uj)duj 

d ., , , , Luji(t,uj)duj 

ot j^ujnQ[uj)auj 

where s{t,uj) + i{t,uj) = nQ{uj) for any t, and no{uj) is a given density which specifies 
probability density function of contact distribution. Using the property that i(t,Lj) = 
no(a') — s{t,uj), we obtain 



d . . . . 

—s{t,w) = -rujs{t,u) 



ujs{t, uj)du} 
Ljno{Lj)du 



(7) 



Models (l2])-(jl!), (E]), and (171) are infinite dimensional dynamical systems. The spe- 
cial form of the models, however, allows us to use well developed tools of the theory of 
heterogeneous populations, which are presented in the following section. 



3 Some facts from the theory of heterogeneous pop- 
ulations 

Here we present some results from the theory of heterogeneous populations in the form 



suitable for our goal noting that more general cases can be analyzed [13|. For the proofs 



we refer to [15[, where similar models are considered. 

Let us assume that there are two interacting populations whose dynamics depend on 
trait values ui and uj2 respectively. The densities are given by ni(t, ui) and n2(t, UJ2), and 
the total population sizes A^i(t) = J^^riiit, u) dui and N2{t) = J^^n2(t,u!) duj2- Obviously, 
more than two populations can be considered, or some populations may be supposed to 
be homogeneous. Assume next that the net reproduction rates of the populations have 
the specific form which is presented below: 

d 

—ni{t,ui) = ni(t,t^i)[/i(vi) + ifi{ui)gi{-vi)], 

f (8) 

—n2{t,UJi) = ra2(t,CJ2)[/2(v2) + 92(^2)], 

where vi = (A^i, A^2, '^2(^)), V2 = (A^i, iV2, <^i(t)), (pi{uji) are given functions, ipi{t) = 
(pi{u!i)pi{t,uji) doji are the mean values of ipiioji)^ and pi(t,uji) = riiit , Ui) / Ni(t) are the 
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corresponding pdfs, i = 1,2. We also assume that ipi{uji), considered as random variables, 
are independent. The system plus the initial conditions 

ni{0,u;i) = N,{0)pi{0,u;i), I = 1,2, (9) 

defines, in general, a complex transformation of densities ni{t,uOi). An effective approach 
to analyze models in the form (El) was suggested in [lH (examples of model analysis are 
given in H Q E, [2H ) . 
Let us denote 

Mi{t,\)= [ e^^'^^^^pi{t,uJi)diUi, 1 = 1,2, 



the moment generating functions (mgfs) of the functions ipi{uJi), Mj(0, A) are the mgfs of 
the initial distributions, i = 1,2, which are given. 

Let us introduce auxiliary variables qi{t) as the solutions of the differential equations 

dqi{t)/dt = gi{yi), g.(0) = 0, 2 = 1,2. (10) 

The following theorem holds 

Theorem 1. Suppose that t G [0,T), where T is the maximal value oft such that (^-(^ 
has a unique solution. Then 

(i) The current means of (fi{uji), i = 1,2, are determined by the formulas 



rfM,(0,A) 



dX 

and satisfy the equations 



^^^^(t) = gi{vi)a^it), (12) 

where afit) are the current variances of ipi{t,uji) , i = 1,2. 

(ii) The current population sizes Ni{t) and N2{t) satisfy the system 

j^N,{t) = N,{t)[f,{^r,) + ^,{t)gi{^r,)], t = 1,2. (13) 

From Theorem 1 follows that the analysis of model ([8])- ([9]) is reduced to analysis of 
ODE system fllOI) . fllip . fll3p . the only thing we need to know is the mgfs of the initial 
distributions. 

Concluding this sections we note that, with obvious notation changes, models (I2l)-(jll) 
and (l5l) fall into the general framework of the master model (IE]). 
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4 Model analysis 



We start with the model 
form 



(jll), which, according to Theorem 1, can be written in the 



Sit) = -mS{t)Iit), S{0) = So, 



dt 

j^m = ms{t)i{t), /(o) = /o, 



d , , 

m 



-I{t), g(0)=0, 
dM{0, A) 



dX 



M(0, A) is a given mgf of ps{0,uj). 

Proposition 1. Model ([Hj) is equivalent to the following model: 

d 



dt 
d 



Sit) = -hiSit))Iit), SiO) = So, 



where 



-lit) = hiSit))Iit), /(O) 



hiS) = So 



'dM~\0,O 




d^ 


5=5/So_ 



and M'^iO,^) is the inverse function to mgf MiO, A). 

Proof. The first equation in ffT^ can be rewritten in the form 



Pit) can be represented as Pit) 



nnM(0,A) 
dX 



\=q{t) 



which gives 



(14) 



(15) 



(16) 



d In Sit) d , , ^ , , , , 

or, using the initial conditions S'(O) = 5*0, g(0) = 0, 

Sit)/So = MiO,qit)), (17) 

which is the first integral to system ( fT^ . Knowledge of a first integral allows to reduce 
the order of the system by one. Since M(0, A) is an absolutely monotone function in the 
case of nonnegative Piui) ^ 0, then it follows that 

qit) = M~'iO,Sit)/So), (18) 
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where M-i(0, M(0, A)) = A for any A. 
Putting (HHD into ([HD gives 

d_g^^^_ dM{0,\) 



SoHt), 

A=A/-i{0,5(t)/So) 



dt d\ 

or, by the inverse function theorem, ( fTSl) with ( fT6l) . □ 



Note that model ([S]) can be reduced to four-dimensional system of ODEs, which, in 
its turn, can be simplified to two-dimensional system. The proof is as in Proposition 1. 
Formally, we have 

Proposition 2. The model ([5]) is equivalent to the model 

pit) = -h,{S)h{I), 

j^m = h,{s)h{i), 

where hi{x), i = 1,2 are given by (fTUj) . 

Combining together Propositions 1 and 2 we obtain the main result of the present 
note. 

Theorem 2. A heterogeneous SI model in the form (^-@, or in the form which 
both describe the contact process with the help of the law of mass action and model het- 
erogeneities in disease parameters such as susceptibility to a disease or infectivity of an 
individual, are equivalent to a homogeneous SI model with a nonlinear transmission func- 
tion. 



An analogous conjecture was made in [27|], where a substantially more complex model 
is analyzed. The strength of Theorem 2 is that it provides an explicit form for the 
nonlinear transmission function. 

Consider a standard gamma distribution with parameters k and u: 

k 

p(0, uj) = -^u'^-^e-"'^, uj ^0, k> 0, u > 0. (19) 
r(A;) 

Let us assume that (3{uj) = uo. The mgf of gamma-distribution is then 

M(0,A) = (1 - A/z/)-^ 
Using Proposition 1 we obtain that 



h{S) = ^ 



s_ 

5^ 



i/fc 



(20) 



From ( l20i) it immediately follows 
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Corollary 1. The power relationship ([T]) with q = l^p = 1 + 1/k can he obtained as a 
consequence of the heterogenous SI model with distributed susceptibility when the initial 
distribution is a gamma- distribution with parameters k and v. 

Corollary 2. The power relationship ([T]) with q = l + l//c2,p = l + can be obtained 
as a consequence of the heterogenous SI model with distributed susceptibility and infectivity 
when the initial susceptibility distribution is a gamma- distribution with parameters ki and 
Vi, and the initial infectivity distribution is a gamma- distribution with parameters k2 and 

Summarizing, we provided a mechanistic derivation of the power law transmission 
function, which was used phenomenologically in many epidemiological models, in the 
case when heterogeneity parameters PiQ exceed one. Originally, these exponents were 
considered to be less than one (e.g., in |23] they are put in the form p=l — a, g = l— 6), 
but no comparison with real world data was provided. 

There is no universal agreement on the values of parameters p, g in ([T]). In ^2J] these 
parameters were estimated when the incidence rate was inferred from epidemic simulations 
on random networks with different degree distributions. In all experiments values of p 



and q were estimated to be less than 1. In contrast to the last observation, in [26|, where 
the transmission function has the form T(S', /) = fiS'^I , it was argued that the exponent p 
should be greater than one. Fitting the solutions of the mean field model with nonlinear 
transmission function into the data obtained from large-scale simulations, it was found 
that p can range from 1.6 to 2. 

In any respect, the question of deriving the power law transmission function on a solid 
mechanistic bases for the case p, g < 1 remains open, whereas the case p, g > 1 is fully 
covered by Corollaries 1 and 2. 



5 Model III and separable mixing 



We rewrite equation ([7]) in the form 



-rujsit^ ijS) 



oj{t)S{t) 
K 



(21) 



where K is the number of contacts, which are made by the total population, ujif) is the 
average number of contacts made by one susceptible individual at time t. We note that 
formally eq. fl2Tl) is not covered by Theorem [T|, because its growth coefficient depends on 
the average parameter value ujit). However, it is possible to extend the theory presented in 
Section 3 to such cases with minor changes in notations (Karev, personal communication). 
In particular, it is possible to show that equation ( l2Ti) is equivalent to the following 
ordinary differential equation: 



d_ 
di 



S{t) = -rh{S) 



1 - 



h{S) 
K 



(22) 
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where h{S) is given by (|TU|) . 

It is interesting to note that we can compare solutions of with solutions of the 
system of ODEs, obtained as a result of large mixing rates in the model on dynamic 
contact network |128|]. For Sl-model system (2.22)-(2.23) from the cited work reads 

= -rMjO, 

where g{x) is the probability generation function for the distribution of the number of 
contacts in the population (this is PGF for pdf rioiu)/ f^no{uj)du!); 9(t) is the fraction 
of individuals that have only one contact and still susceptible by the time t; r is the 
transmission rate; and Mj is the fraction of contacts made by infected individuals. The 
number of susceptible individuals is given by S(t) = g{6{t)). 

To compare models (l22i) and (l23ll we need to specify the initial conditions. Since 
model (l23l) deals with PGF of the number of contacts of the total population, and eq. 
fl22|) incorporates mgf of the number of contacts of susceptible individuals it is reasonable 
to expect some discrepancy of the corresponding solutions if we use the same pdf for these 
purposes. See Fig. [T]for three solutions. 

As can be seen from Fig. [T]the best agreement os found when we use A = 1, i.e., the 
average number of contacts equals to 1. In this case two solutions coincide. In the cases 
A > 1 or A < 1 there is some divergence, although the limiting behavior of the models is 
the same. 
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1000 




1 2 3 4 5 

t 

Figure 1: Comparison of the solutions of system (l23ll (red line) with solution of eq. ( |22l) 
(yellow bold line). Poisson distribution was used with parameters A = 1, 1.5, 0.5 from 
top to bottom, r = 2. It was assumed that the population size is = 1000. The 
initial conditions for problem ( !23l) were chosen such that 6{0) = 1 — e, M/(0) = e, where 
e = 0.01. S{0) for ( l22l) was found as g{l — e). The dotted line shows Npo, where po is 
the proportion of individuals in the population who do not make the contacts 
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